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^D Abstract. Plotting spectra of a range of almost Mathieu operators reveals a 

"^ beautiful fractal-like image that contains multiple copies of a butterfly image. We 

i^ demonstrate that plotting the butterflies using a gap-labelling scheme based on K- 

, theory or Chern numbers reveals systematic discontinuities in the gap positioning. 

|I^ A proper image is produced only when we take into account these discontinuities, 

*^ and close the butterfly wingtips at the points of discontinuity. A conjecture is 

^H presented showing a simple formula for locating the discontinuities, and numerical 

f — evidence is given to support the conjecture. We also present new renderings of this 

butterfly. 

< 

q 

^ 1. Introduction 






The inspiration for this work is the diagram in Figure 1, a fractal-hke structure 



known as the Hofstadter butterfly (Hofstadter, 1976), which represents the energy 



levels of an electron travelling through a periodic lattice under the influence of a 



magnetic field; that is, a Bloch electron (Brown, 1964). Each horizontal line in the 



^H image corresponds to union of spectra of certain bounded self-adjoint operators Hgif, 

^ on the Hilbert space /^(Z), defined by the equation 

^ (-f^e,</,x)„, = x„_i + x„+i + 2 cos(27m6' + 0)x„, 

(^ where the parameter 9 G [0, 1] corresponds to the vertical position in the image 

O in Figure 1, and is a physical phase parameter. This operator, known as the 



almost Mathieu operator and also as Harper's operator (Last, 1994), is a discrete, 
.^ one-dimensional Schrodinger operator, which has been the object of study for many 

^ important physical problems (Goldman, 2009). 



j^ Equivalently, each horizontal line in Figure 1 represents the spectrum of the self- 

adjoint operator 

Hq = U + U* + V + V* 

where unitaries u, v are generators of the rotation C*-algebra Ae, characterized by a 



commutation rule 



vu = e^^'\v. 
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Figure 1. The Hofstadter butterfly. 



These particular algebras are known as non-commutative tori and are of considerable 



interest in noncommutative geometry (Connes, 1994). The Cantor-like structure of 
the spectra for operators in these algebras has been a topic of study in C*-algebras 



Simon, 1982; Choi et al., 1990; Last, 1994 Puig, 2004) 



for quite some time, as discussed in (Avila and Jitomirskaya, 2006 Bellissard and 
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Figure 2. Drawing only the butterflies. 



Our interest in this paper is in the mechanics of drawing the image in Figure 1, its 



related image in Figure 2, and what this reveals about gap labelling theorems (Bel- 



lissard, 1990 Goldman, 2009 Kaminker and Putnam, 2003 Ypma, 2007) for these 



Schrodinger operators. This particular image in Figure 1 was created in high reso- 



lution, using the numerical package MATLAB (TheMathWorks, 2008) to accurately 
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compute the eigenvalues of certain tridiagonal matrices, then plotting line segments 
between eigenvalues using PostScript commands to ensure a high quality image. It is 



useful to refer to (Casselman, 2005) for details on how to create such high resolution 
mathematical illustrations in PostScript. 

In the finely rendered image in Figure 1, the viewer's eye is immediately drawn to 
butterfly-like structures in the diagram which repeat over and over at various scales 
and positions. Following Casselman's dictum that mathematical drawings should 



represent concepts, not accidents (Casselman, 2005), it is fascinating to consider how 



we might draw only the butterflies. In Figure 2, we do exactly that, drawing precisely 
the butterflies that are inherent in the image in Figure 1, without drawing all the 
spectral lines themselves. 

This new high resolution image was also created using Postscript commands, based 
on numerical calculations of spectra of operators generated in MATLAB. This new 
image highlights the fractal symmetry of the diagram and suggests further investi- 
gation on the explicit symmetries. 

The goal of this paper is to record the algorithmic process used to image only 
the butterflies, using the gap labelling procedure for the almost Mathieu operator. 
We discover through this process that, from a certain point of view, there are dis- 
continuities in gap labelling, and we present a conjecture for specifying where these 
discontinuities occur. It appears this conjecture is new to the literature, but we are 
at a loss to prove it. The conjecture is verified for low rational values of parameter 
9. 

A secondary goal is to reveal this new approach to rendering a high resolution 
image of the Hofstadter butterfly. 

2. Spectral calculations 
The spectral calculations for kg are well-known, and are particularly straight- 



forward when 9 is rational (Choi et al., 1990 Lamoureux and Mingo, 2007). For 



9 = p/q in reduced form, irreducible representations of the rotation C*-algebra Ag 
are obtained using q x q matrices U, V, with U a cyclic permutation matrix and V 
a diagonal matrix with entries given as consecutive integer powers of e^'^*^'"'. Under 
the irreducible representation, the operator he maps to the q x q matrix 

H0,,,,z, =ziU + ziU* + Z2V + ziV\ 

where Zi, Z2 are complex numbers of modulus one, specifying the different irreducible 
representations. The spectrum of hg is just the union of the point spectra for the 
various matrices Hg^zi,z2 over the range of zi, Z2 

The characteristic polynomial of the qxq matrix Hg^zi,z2 is invariant under changes 
of parameters Zi,Z2, except for a zero-order term zf + z^'^ + Z2 + Z2^ . Consequently, 
the union of spectra over the set of parameters zi, Z2 is just the inverse image, under 
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a degree q polynomial, of the range of values for the term zf + Zi'^ + Z2 + -2^^. This 
range is the interval [—4, 4] , and the inverse image is a set of q distinct intervals on 
the real line. 

These q intervals are all disjoint, except in the case q even, in which case the 
central two intervals just touch at the spectral point zero. Thus, there are precisely 
g — 1 gaps between intervals in the spectrum of hg, where we include the "empty 
gap" at spectral value zero, for convenience. 

To compute the spectral intervals, it suffices to compute the endpoints of each 
interval and then connect with a line. The endpoints of the interval are obtained by 
choosing values of zi, Z2 such that the constant term zf + z^'^ + -z^ + ^2'^ obtains its 
extreme values ±4. 

The q X q matrix i/e,^i,z2 is almost tridiagonal - and by choosing the extreme 
values for -21,2:2, we can apply explicit Givens rotations to obtain a preferential tri- 



diagonal form for Hg^zi,z2i as described in (Lamoureux, 1997). In fact, the matrix then 



decomposes as the direct sum of two tridiagonal matrices, each of size approximately 
q/2 X q/2. Finding eigenvalues for these tridiagonal matrices is an order g^ operation, 
which is a fast computation in this application. 

3. Gap labelling: inverse slope 1,-1 

To draw the butterfly wings, we do not draw the spectral lines, but instead draw 
curves around the spectral gaps. 

The spectral gaps are conveniently labelled using a Diophantine equation, where 
integer parameters t, s are fixed, and we select the r-th gap in spectrum 6 = p/q 
using the formula 

r = t*p — s*q. 



This indexing is known as gap labelling, as described in (Bellissard, 1990 Goldman 
2009 Kaminker and Putnam, 2003; |Ypma, 2007), and the parameters (t, s) are 



related to indices in K-theory, or equivalently, Chern numbers. Here, the point is 
that at each spectral level 9 = p/q, the r-th gap does actually appear (if we also 
include the zero-length gap at spectral value 0), and we have a convenient way of 
marching through all the gaps. 

Of course one must choose reasonable values oi t,s in order to give a reasonable 
r value in the range of 1 ... g — 1, since there are only q — 1 gaps. Thus for integers 
t > we require 

<s < t- 1 

and for t < we require 

t < s < -1. 
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Note we are using a slightly different convention from the literature for our gap 
labelling, introducing a minus sign into the Diophantine equation 

r = t*p — s*q. 

This allows us to keep indices (t, s) with the same sign, which is just a matter of 
convenience in our algorithm. 

As an example, we set t = 1 which forces s = 0, and we plot the gaps in first image 
shown in Figure 3. There are two curves in this first image. First, the curve on the 
left of the butterfly wing, which consists of the right endpoints of spectral intervals 
to the left of the gap. Second, the curve on the right of the butterfly, consisting of 
the left endpoints of the spectral intervals to the right of the gap. 

For the negative slope t = — 1, we are forced to take s = — 1, and we see the plot 
for the downward sloping butterfly wings, also shown in Figure 3. 





Figure 3. Butterfly wings at slope 1 and -1. 

It is useful to refer to the parameter t as the inverse slope, for the Diophantine 
equation r = t*p — s*q becomes, after division by q, a linear equation 



t* 



p 



s, or X 



t*9-s, 
q q 

where t indicates the slope between variables x and 9. 

Notice something special happens at the wingtips, where 6 = 0,1. For these values 
of 9, there is no spectral gap, as the spectrum consists of the single interval [—4,4]. 
So it does not make sense to identify the r-th gap there. From Figure 1, though, it 
is clear what should happen: the butterfly wing should close to a tip. To get this 
correct, we must include the endpoints of the interval [—4, 4] as the closing point for 
the wing. 

We will see in the next sections that we need additional closures to complete the 
smaller butterfly wingtips. 
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4. Gap labelling: inverse slope 2 

We continue with the gap labeUing, looking at gaps with label r = 2*p — s*q\ 
that is, wings with inverse slope 2. For s = 0, 1, we get the wings shown in Figure 4, 
the first diagram covering a 6 range of [0, 1/2] and the second covering range [1/2, 1]. 
Note that we must close the wingtips at 6* = 0, 1/2, 1, just as we did in the last 
section, selecting the min/max spectral endpoints, rather than a gap. 

However, we also see a new phenomena. Again in Figure 4, on the left diagram, 
there is an extraneous segment connecting two distinct wings, at the vertical value 
9 = 1/3. Similarly, in the the right diagram, there is a segment at ^ = 2/3. This is 
in fact a discontinuity in the gap labelling. At ^ = 1/3, there is a jump in the gap. If 
we try to draw that jump line, we would get a hne cutting across the larger butterfly 
wing drawn in the last section. (The wing with inverse slope t = —1.) 

So, we should remove these discontinuities. Again, we want to close the wingtips, 
so at ^ = 1/3 we choose the same spectral endpoint to close the wingtip. Specifically, 
when coming up from below, we pick a right endpoint of a spectral interval. When 
coming down from above, we pick a left spectral endpoint. This is exactly as we did 
with the endpoints ^ = 0, 1. 

With these discontinuities removed, we see proper wingtips, shown in Figure 5. 





-3 -2 



Figure 4. Butterfly wings at inverse slope 2, s = 0, 1. Note the 
cutting segments at 6* = 1/3, 2/3. 



5. Gap labelling: inverse slope 3 

We continue the construction of butterfly wings, by labelling gaps with inverse 
slope t = 3. In this case we have three choices for s, as s = 0, 1, 2, giving the 
diagrams in Figures 6, 7, 8. Again, Figure 6 only covers the range of 6 in [0, 1/3], 
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Figure 5. Butterfly wings at inverse slope 2, s = 0, 1. Cutting seg- 
ments removed. 

Figure 7 covers the range [1/3,2/3] and Figure 8 the range [2/3,1]. We must close 
wingtips at the endpoints of these ranges (i.e. at 0, 1/3,2/3, 1), but we also notice 
cutting segments again. 





Figure 6. Butterfly wings at inverse slope 3, s = 0. Cutting segments 
on the left, at 6 = 1/5, 1/4, and removed on the right. 

In Figure 6, on the left, the cutting segments are clearly at 6' = 1/5 and 6 = 1/4. 
Removing the cutting segments and closing the wingtips gives the right diagram in 
Figure 6. 

In Figure 7, on the left, the cutting segments are not quite so clear, but a close 
examination shows there are segments at ^ = 2/5 and 6 = 3/5. Removing the cutting 
segments and closing the wingtips gives the right diagram in Figure 7. 
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Figure 7. Butterfly wings at inverse slope 3, s = 1. Cutting segments 
on the left, at 9 = 2/5, 3/5, and removed on the right. 





Figure 8. Butterfly wings at inverse slope 3, s = 2. Cutting segments 
on the left, at ^ = 3/4,4/5, and removed on the right. 

In Figure 8, on the left, the cutting segments are at 6' = 3/4 and 6 = 4/5. Removing 
the cutting segments and closing the wingtips gives the right diagram in Figure 8. 



6. Gap labelling: inverse slope 4 

As we go to higher inverse slope t, it is harder to eyeball where the discontinuities 
should be. 

For some ranges of 6, the pattern is clear. For t = 4, there will be four ranges of 
6, namely 

[0,1/4], [1/4,1/2], [1/2,3/4], [3/4,1]. 
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In the first interval, we expect wingtips (or discontinuities) at 

e = 0, 1/7, 1/6, 1/5, 1/4. 
In the last interval, we expect wingtips (discontinuities) at 

= 3/4,4/5,5/6,6/7,1, 
which we can guess just from the patterns that occurred at slopes t = 2, 3. 
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X; 0.1403 
Y: 0.4289 

■ " 




X: 1.539 
Y: 0.4005 








■ 






X: 1,306 
Y: 0.3341 








1 






■ X: 0.02254 
Y: 0.2861 

■ 







Figure 9. A plot of jumps at inverse slope 4, big ones labeled. 



However, on the interval [1/4,1/2], it is not so clear, so it helps to do some 
numerical experimentation. In Figure 9, we plot the jumps in the gap as 6 in- 
creases. At the labelled points, we note significant jumps at the vertical values 
9 = 0.2861, 0.3341, 0.4005, 0.4289. 

Estimating in terms of fractional values, these jumps are at 6 values 2/7, 1/3, 2/5, 3/7. 
Significantly, there are four jump values - whereas from the other examples, we might 
have only expected 3 = t — 1 values. In Section 8, on the gap labelling conjecture, 
we show why there are four jump points. 

Using this jump data, we remove cutting lines and obtain the inverse slope 4 wings 
shown in Figure 10. 
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Figure 10. Butterfly wings at inverse slope 4, s = 0, 1, 2, 3. 

7. Sanity check: another jump at inverse slope 4? 

Looking carefully at Figure 9, it appears there is another signiflcant gap jump at a 
9 value of about .38. Examining our numbers carefully in MATLAB, we flnd a value 
of ^ = 0.3747, which is about 9 = 3/8. 

Is there a jump here? No. This value corresponds to the spectral point zero, the 
centre of the butterfly pattern, where the inverse of the butterfly curves flattens out 
- so we have large movement, even though there is no jump. 

We note that t*9 — s, ai 9 = 3/8, gives the value 1/2, the midpoint of the butterfly, 
or spectral value of zero. 



8. Gap labelling: inverse slope 5 

At inverse slope 5, it is very difficult to guestimate where the discontinuities should 
be. However, in the next section of this paper, we have a conjecture predicting the 
gap jumps for 9 values of the form 

S + So 



9 



t + to 



, for any < i < to, < s < t. 



where we are following the gap curve r = to*p — So*q, with to > Sq > 0. For inverse 
slope 5, we set to = 5 and So = 0, 1, 2, 3, 4. We use these fixed values, and the range 
of t, s allowed above to pick out discontinuities, and then plot. Again, we must close 
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the wingtips following the usual algorithm of matching endpoints of the appropriate 
spectral intervals. 

Figure 12 shows the collection of all the butterfly wings with inverse slope 5, all 
cutting lines removed, and wingtips closed, using the conjectured gap discontinuities. 
The image looks correct, suggesting that the conjecture works for inverse slope 5. 



.^ 



/ 



,^ 



^ 



.^ 




For r = rp + s'q, 5<= t <= 5 



2 3 



Figure 11. Butterfly wings at inverse slope 5, s = 0, 1, 2, 3, 4. 



9. Gap labelling: discontinuity conjecture 

We observe the apparent discontinuities in the gap labelling can be determined by 
the following algorithm. 

Algorithm: For a given labelling r = to * p — SoQ, say with to > 0, we plot the 
positive sloping line in [x, 9) given by the equation 

X = to* 6 — Sq, 

and all the negative sloped lines in (x, 9) given by 

X = —t* 9 + s, 

where (t, s) are integers satisfying < s < t < to. Note that (— t, —s) is a valid index 
for a negative slope gap labelling. 

Conjecture: The intersection of all these lines, in the region < x < 1 give the 
gap discontinuities. 
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Visually, we can see the intersections. For instance, in the case of a line slope 
to = 2, we see in Figure 12 the relevant lines for Sq = (on the left) and Sq = 1 on 
the right, giving the intersections at ^ = 1/3 and 6 = 2/3 respectively. 





Figure 12. Line intersections specify the discontinuity points for 
slope 2, s = 0, 1. On left, red line indexed as (2,0) has intersection at 
6 = 1/3. On right, red line indexed as (2,1) has intersection at ^ = 2/3. 

In the case of a line slope to = 3, we see in Figure 13 the relevant lines for 
So = 0, So = 1, and So = 2 has the intersections at 6* = 1/5, 1/4, 9 = 2/5, 1/2,3/5, 
6 = 3/4,4/5, respectively. The intersection at 6 = 1/2 is not really a discontinuity, 
as the spectral intervals close here (spectral point zero). However, it does not hurt 
the algorithm to include this point as a jump point. 






Figure 13. Line intersections specify the discontinuity points for 
slope 3, s = 0, 1, 2. Left, red line at (3,0) has intersections at 6* = 
1/5, 1/4. Centre, red line at (3,1) has intersections at ^ = 2/5, 1/2, 3/5. 
Right, red line at (3,2) has intersections at 6* = 3/4, 4/5. 
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Figure 14. Line intersections specify the discontinuity points. Red 
line indexed as (4,1) has four intersections, at 6* = 2/7, 1/3, 2/5, 3/7. 

The case of slope to = 4 is particularly interesting, as this is the first level where 
the "pattern" of jump discontinuities is not immediately clear. Looking at Fig- 
ure 14, for the line of index (4,1), we see gap discontinuities predicted at values 
6 = 2/7, 1/3, 2/5, 3/7, just as we determined numerically in Section 6. This seems 
to confirm that the algorithm is promising (and perhaps correct). 

It is useful to state specifically the conjecture about the gap discontinuities: 

Theorem 1 (Conjecture). For gaps labelled by the Diophantine equation 

r = to* p — So* q, with to > So> 0, 
there are discontinuities in the gap labelling at values 9 of the form 

Q^S_o±^ 
to+t' 

for all integers t, s satisfying < s <t < to, and restricted to the open interval 

So + S fso So+l 
u = t 



Closing the gaps at the discontinuities gives the butterfly wingtips. 

We do not have a proof of this result, only numerical evidence to confirm it. The 
predicted discontinuities also include the zero-length gaps in the spectral intervals. 
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where two intervls touch at spectral value zero. Including these points in the plotting 
algorithm does not disturb our images. 

As an example of implementing this theorem, with {to,So) = (4,1), we expect 

discontinuities at 

1 + 1 1 + 1 1 + 1 1 + 2 

4 + 1' 4 + 2' 4 + 3' 4 + 3 
which gives exactly the points 9 = 2/5, 1/3, 2/7, 3/7 as expected. 

As the final numerical test of the conjecture, we note that the butterfly image 
in Figure 2 was computed using this algorithm for detecting discontinuities in the 
gap labelling, and we observe that Figure 2 is an excellent rendering of the butter- 
flies apparent in the original spectral map in Figure 1. Figure 2 includes all gap 
labelling with all non-zero inverse slopes in the range —10 < t < 10. It looks correct, 
suggesting the algorithm to remove gap discontinuities is working correctly. 

In the appendix we include a flgure with gap labelling for all non-zero inverse 
slopes in the range —20 < t < 20. Again, we see no obvious breakdown in the 
algorithm used to remove the gap discontinuities. 

10. Sanity check: zooming in on discontinuities 

As a check that our algorithm for predicting the discontinuities in the gap labelling 
is correct, we will examine closely the image at an unexpected discontinuity. 

As mentioned in Section 6, the discontinuity on line (to, Sq) = (4, 1) at ^ = 2/7 was 
unexpected. For this line, we were expecting only three discontinuities, as there are 
big jumps at the three values 9 = 1/3, 2/5, 3/7. The value 2/7 doesn't quite flt the 
apparent pattern, and yet the conjecture says there should be a discontinuity there. 

To test this, we plot side-by-side the butterfly image using spectral lines, and 
using the new algorithm for the drawing the labelled gaps, zoomed in at the possible 
discontinuity point 9 = 2/7. The two images are shown in Figure 15. Notice how 
the horizontal line exactly at the center (left image) has a gap - so there really is a 
spectral gap at ^ = 2/7. The image on the right, showing the butterfly, accurately 
captures this gap - that is, we see the main butterfly image here has a gap at the 
centre. So the algorithm successfully identifled the gap discontinuity in this case. 

It is signiflcant, and worth noting, that the two wingtips at the centre of the dia- 
gram do not meet, and are slightly offset. It is an interesting phenomenon, accurately 
rendered both with spectral lines (left image in Figure 15) and with the butterflies 
(right image in Figure 15). 

11. Sanity check: other coupling constants 

It is a useful check on the algorithm to verify that the butterfly structures maintain 
the same basic form under small modiflcations of the setup of the problem. The Bloch 
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Figure 15. Comparing the discontinuity at 6 = 2/7, zoomed in. 



electron model includes a physical parameter called the coupling constant A which 
is generically set equal to 2. For the general problem, one can consider self-adjoint 
operators of the form 

he = u + u* + -{v + V*), 

where A is a positive real parameter. It is known that the gap structure for the 



corresponding spectra remains the same (Bellissard and Simon, 1982). 



We run our code as before, only changing the coupling parameter A and see that 
the basic form of the result looks the same. Figure 7 shows three results, for three 
different values of the coupling parameter. Note the width of the result varies with 
A, and in fact the width is proportional to 2 + A. Other than the width differences, 
the butterfly structures look very similar, and we do not observe any breakdown in 
the algorithm that predicts the gap discontinuities. Again this is evidence of a good 
conjecture. 



12. Future work 

It would be useful to have a function that can select, then draw, individual pairs 
of butterfly wings. As things stand now, we only get half-wings of certain slopes, 
and it is not clear how to label them and put them all together. 

This may be related to the fractal structure of the Hofstadter butterfly, which we 
supsect is somehow indexed by 5^2 (Z). 
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Figure 16. Three versions with couphng constant A = 1.6,2.0,3.0, respectively. 

13. Summary 

We have presented a new rendering of the Hofstadter butterfly based on explicit 
drawings of the spectral gaps using the gap labelling method of K-theory. Disconti- 
nuities in the gap positions are revealed by this method, and we present a conjecture 
for identifying where the discontinuities occur, and how to "close" the the wing tips 
at the discontinuities. We have presented numerical evidence to suggest the conjec- 
ture is correct, and results in an accurate rendering of the fract-like structure for the 
Hofstadter butterfly. 

As a side note, we observe the utility of programming our mathematical illustra- 
tions directly in the PostScript language, for optimum resolution. 
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Appendix 1: Code 

We list the MATLAB code that draws the butterfly wing. 

There is a call to a function Heigs{p, q, 2) which produces the list of eigenvalues 
for the spectral lines of the almost Mathieu operator. This function is discussed 
elsewhere. 

7o script diophcLntS.m 

7o MATLAB code to draw a range of butterfly wings 

7. 

7o We will draw a whole rgmge of butterfly wings, from tmin to tmax 

7o Follow the diophantine parameterization for the gap count: 
7o r = t*p - s*q. 

7o We run through a range of t values, 1 <= t <= 5, say. 
7. And <= s < t 

tmin = 1 ; tmax = 5 ; 

q_max = 50; 7o maximum denominator we want to consider 

elf; hold on 7 clear the figure and set hold on 

for t=tmin:tmax 
for s=0: (t-1) 

7o first we compute the jump points 
start_pt = s/t; 
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end_pt = (s+l)/t; 
jmp_pts = [start_pt, eiid_pt] ; 
for tl=l:(t-l) 
for sl=l:tl 

new_pt = (sH-s)/(tl+t) ; 

if (new_pt > start_pt && iiew_pt < end_pt) 

if (sum(new_pt == jmp_pts)==0) 7, not in list yet 

jmp_pts = [jmp_pts new_pt] ; 
end 
end 
end 
end 

jmp_pts = sort (jmp_pts) ; 

°/o Now we do the butterfly wings, from lower left to upper right 

for k = 1 : (length(jmp_pts)-l) 

start_pt = jmp_pts(k); "/ set up the start and end points, plot between 
end_pt =jmp_pts(k+l) ; 

theta_list = [] ; 7, the list of various points to plot 
lef t_list = [] ; 7o the left side of the butterfly gaps 
right_list = [] ; 7 the right side of the butterfly gaps 

for q=l:q_max 

for p=0:q 7 we use the fact that gcd(0,l) = gcd(l,l) = 1 to get ends 
r = t*p - s*q ; 

if (gcd(p,q)==l && start_pt <= p/q && p/q <= end_pt ) 
if (start_pt < p/q && p/q < end_pt) 
X = Heigs(p,q,2) ; 

theta_list = [theta_list ,p/q] ; 7. add in a new theta value 
left_list = [left_list,x(2*r)] ; 
right_list = [right_list,x(2*r+l)] ; 
end 
if (p/q==start_pt) 

X = Heigs(p,q,2) ; 

theta_list = [theta_list ,p/q] ; 7. add in a new theta value 
left_list = [left_list,x(2*r+l)] ; 
right_list = [right_list,x(2*r+l)] ; 
end 
if (p/q==end_pt) 

X = Heigs(p,q,2) ; 

theta_list = [theta_list ,p/q] ; 7 add in a new theta value 

left_list = [left_list,x(2*r)] ; 

right_list = [right_list,x(2*r)] ; 
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end 

end 
end 
end 

[theta_list, IX] = sort (theta_list) ; 
left_list = left_list(IX) ; 
right_list = right_list(IX) ; 
plot(left_list,theta_list, '-') 
plot (right_list , theta_list , ' - ' ) 
plot (-lef t_list , theta_list , ' - ' ) 
plot (-right_list , theta_list , ' - ' ) 

xlabel(['For r = t*p - s*q, ' ,num2str (tmin) , '<= t <= ' , num2str(tmax)] ) 
xlim([-4,4]);ylim([0,l]); 

end 7o we end the k loop 
end 7o we end the s loop 
end % we end the t loop 

Appendix 2: EPS file 

Here is a sample of the Encapsulated Postscript file (EPS) used to create the 
butterflies in the high resolution Postscript files. There is a header including the 
format information, then a series of Postscript commands that set up the scaling, 
line width, and finally a series of moveto/lineto command that draws the lines. 

We use MATLAB to create the numbers, then save the list of commands directly 
to a file that can be used by any Postscript reader. 

7. ! PS-Adobe-3 . EPSF-3 . 

7o7oCreator: Michael P. Lamoureux, Copyright 2010 

y.y.Title: Almost Mathieu Butterflies 

y.y.CreationDate: 16-Feb-2010 

7.y.DocumentData: Clean7Bit 

y.y.Origin: 

y.y.BoundingBox : -400 400 800 

100 100 scale 

.0005 setlinewidth 

1 setlinecap 

newpath 

-4.000000000000e+00 .OOOOOOOOOOOOe+00 moveto 

-3.876300213013e+00 1 .600000000000e-01 lineto 

-3.873816420419e+00 1 .632653061224e-01 lineto 

3.625495844745e+00 7.836734693878e+00 lineto 
3.632788171307e+00 7.840000000000e+00 lineto 
4.000000000000e+00 8. OOOOOOOOOOOOe+00 lineto 
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stroke 

newpath 

4.000000000000e+00 0. OOOOOOOOOOOOe+00 moveto 

3.876300213013e+00 1 . 600000000000e-01 lineto 

3.873816420419e+00 1 . 632653061224e-01 lineto 

-3.625495844745e+00 7.836734693878e+00 lineto 
-3.632788171307e+00 7.840000000000e+00 lineto 
-4. OOOOOOOOOOOOe+00 8. OOOOOOOOOOOOe+00 lineto 
stroke 



Appendix 3: more plots 

Now that we have the technology, let us plot two more images for looks. 

The first image, Figure 17 shows 5 levels of butterflies, for a clean, simple im- 
age. The second image. Figure 18, shows 20 levels of butterflies, for a denser, more 
complex image. 

In all cases, the algorithm for predicting the gap discontinuities successfully re- 
moves any extraneous lines. 
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Figure 17. Five levels of butterflies, with denominator q in the range 1...50. 
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Figure 18. Twenty levels of butterflies, with denominator q in the 
range 1...50. 
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